MBD Estadística - Practica 1

A continuación se expone el código y comentarios resultados de la realización de la práctica 1.

0. Imports

Importar librerias

library(tidyverse)
library(tidyr)
library(dplyr)
library(knitr)
library(ggplot2)
library(stringr)
library(reshape)
library(minerva)
library(heatmaply)
library(kableExtra)
library(factoextra)
library(ggbiplot)
library(randomForest)
library(ROCR)
library(ggcorrplot)
library(ModelMetrics)

Importar datos

Aunque se nos entrega el fichero de train entero y podriamos obtener la variable respuesta, se trataran cómo si de un caso real se tratara, y todo el EDA y el training se usara solo el p1_train. El p1_test se usara únicamente para validar los resultados al final de todo de la practica.

train <- read.csv("./data/train.csv", head = TRUE, stringsAsFactors=TRUE)

p1_train <- read.csv("./data/p1_train.csv", head = TRUE, sep=";", stringsAsFactors=TRUE)

p1_test <- read.csv("./data/p1_test.csv", sep=";", head = TRUE, stringsAsFactors=TRUE)

1. Preprocesado de los datos

En este apartado se va hacer el análisi exploratorio de los datos.

head(p1_train)

Eliminamos la variable id del dataset ya que es el indice, y no aporta información.

p1_train <- p1_train %>% select(-id)

Summary de los datos

summary(p1_train)

      year           hour           season         holiday        workingday        weather          temp           atemp      
 Min.   :2011   Min.   : 0.00   Min.   :1.000   Min.   :0.000   Min.   :0.0000   Min.   :1.00   Min.   : 0.82   Min.   : 0.76  
 1st Qu.:2011   1st Qu.: 6.00   1st Qu.:2.000   1st Qu.:0.000   1st Qu.:0.0000   1st Qu.:1.00   1st Qu.:13.94   1st Qu.:16.66  
 Median :2011   Median :12.00   Median :3.000   Median :0.000   Median :1.0000   Median :1.00   Median :20.50   Median :24.24  
 Mean   :2011   Mean   :11.57   Mean   :2.506   Mean   :0.029   Mean   :0.6773   Mean   :1.41   Mean   :20.27   Mean   :23.70  
 3rd Qu.:2012   3rd Qu.:18.00   3rd Qu.:4.000   3rd Qu.:0.000   3rd Qu.:1.0000   3rd Qu.:2.00   3rd Qu.:26.24   3rd Qu.:31.06  
 Max.   :2012   Max.   :23.00   Max.   :4.000   Max.   :1.000   Max.   :1.0000   Max.   :3.00   Max.   :41.00   Max.   :45.45  
    humidity        windspeed          count      
 Min.   :  0.00   Min.   : 0.000   Min.   :  1.0  
 1st Qu.: 46.00   1st Qu.: 7.002   1st Qu.: 41.0  
 Median : 62.00   Median :12.998   Median :145.0  
 Mean   : 61.77   Mean   :12.802   Mean   :191.4  
 3rd Qu.: 77.00   3rd Qu.:16.998   3rd Qu.:283.0  
 Max.   :100.00   Max.   :56.997   Max.   :977.0  

1.0 Nan handling

Se puede apreciar que no encontramos ningún nan en las columnas del dataset, haciendo que no tengamos que tratar con ellos.

p1_train %>% group_by() %>% summarise_all(funs(sum(is.na(.))))

1.1 Duplicates handling

Se eliminan las filas duplicadas del dataset.

p1_train <- p1_train %>%
                distinct()
p1_train

1.2. Separar datos entre numericos y categoricos

Para las variables que son categoricas, aunque su representación sea numerica debemos tratarlas cómo categoricas a nivel dato para obtener el resultado deseado.

Las variables categoricas del dataset segun el enunciado son: workingday, holiday, weather y season.

p1_train <- p1_train %>% mutate(workingday = as.factor(workingday),
                                   holiday = as.factor(holiday),
                                   weather = as.factor(weather),
                                   season = as.factor(season))

2. Exploratory Data Analysis

Eliminamos la variable target del dataset para hacer un estudio de las columnas predictivas. Posteriormente se hará el estudio de la variable target independientemente.

predictors <- p1_train %>% select(-count)

2.0. Hacer boxplots de los campos y mirar outliers en los datos

(Sobre variables numericas)

meltPred <- predictors %>% select(where(is.numeric)) %>% melt()

Using  as id variables

meltPred %>%
 ggplot(aes(factor(variable), value)) +
   geom_violin(width=1, color = "gray", alpha = 0.2, fill = 'green' ) +
    geom_boxplot(width=0.3, color="black", fill='blue', alpha=0.3, outlier.colour="red", outlier.shape=8,
             outlier.size=1, notch=TRUE) + facet_wrap(~variable, scale="free")

Notch went outside hinges
ℹ Do you want `notch = FALSE`?

Estos gráficos representan un boxplot y un violinplot sobre todas las variables numericas del dataset. Podemos apreciar:

• La columna year contiene únicamente dos valores 2011 y 2012.

• La columna hour se ditribuye entre el 0 y el 24, con una media situada en los 12. No encontramos outliers.

• La columna temp y la columna atemp siguen distribuciones muy similares, la media esta alrededor del 25. No encontramos outliers.

• La columna humidity se distribuye entre el 0 y el 100 con una media alrededor de los 60. No encontramos otuliers.

• La columna windspeed se distribuye entre el 0 y el 60 con una media alrededor de los 15. Encontramos numersos outliers cuando el valores es mayor de 30.

2.1 Matriz de correlaciones de pearson de las variables numericas

corrp <- p1_train %>% select(where(is.numeric)) %>%  cor(method = 'pearson')
corrp %>% ggcorrplot(hc.order = TRUE, type = 'lower', lab=TRUE,
                      outline.col = "white",
                       ggtheme = ggplot2::theme_gray,
                       colors = c("#6D9EC1", "white", "#E46726")
                      )

Se puede apreciar cómo la variable atemp y temp estan muy correlacionadas, de ambas, seleccionaremos solamente la columna temp y eliminaremos la columna atemp para una mejora en la prediccion. Con la variable target count, estan directamente relacionadas las columnas temp, year, hour y windspeed y la columna humidity esta indirectamente realcionada.

Eliminamos la columna atemp del dataset.

p1_train <- p1_train %>% select(-atemp)
p1_test <- p1_test %>% select(-atemp)
predictors <- predictors %>% select(-atemp)

2.2 matriz de correlaciones lineales y no lineales con la variable target

Calculamos la correlación no lineal MIC de las columnas con la variable target.

numpred <- predictors %>% select(where(is.numeric)) 
corrs <- data.frame(MIC = mine(numpred, y = p1_train$count, alpha = 0.7)) %>% select(MIC.Y)
corrs$Pearson <- cor(numpred, y = p1_train$count, method = 'pearson')
rownames(corrs) <- rownames(corrs$Pearson)
corrs

Podremos apreciar que las correlaciones no lineadas calculadas son (en proporcion) similares a las correlaciones de Pearson.

cbind(corrs$MIC.Y, corrs$MIC.Y ) %>% heatmaply_cor(xlab = "",
              ylab = "Columnas", main = "Correlacion MIC de las columnas numericas con la variable target 'count'", cellnote = cbind(corrs$MIC.Y, corrs$MIC.Y ), k_col = 1, k_row = 1)

NA

El gràfico superior tiene 2 columnas pero se debe considerar una sola (estan repetidas), ha sido la unica manera que hemos encontrado de hacer el gráfico del MIC.

2.3 PCA

Realizamos un PCA sobre los datos para extraer los componentes principales y estudiarlos.

pca <- predictors %>% select(where(is.numeric)) %>% sample_n(200) %>%
  prcomp(scale. = TRUE, center = TRUE)
summary(pca)

Importance of components:
                          PC1    PC2    PC3    PC4    PC5
Standard deviation     1.2906 1.0555 0.9785 0.8406 0.7458
Proportion of Variance 0.3331 0.2228 0.1915 0.1413 0.1112
Cumulative Proportion  0.3331 0.5560 0.7474 0.8888 1.0000

pca %>% 
  ggbiplot::ggbiplot(scale = 1)

Gràfico que representa la ponderación de las variables utilizadas para el PCA.

3. Analisis exploratorio de la variable target ‘count’

La variable target ‘count’ es una variable numerica que representa el numero de bicicletas alquiladas en las grandes ciudades. La variable tiene una media de 190 y un std de 182, un valor bastante elevado que indica mucha variación en la variable.

corrs$Pearson

                [,1]
year       0.2623084
hour       0.4037404
temp       0.3971578
humidity  -0.3245222
windspeed  0.1078587

Correlaciones de pearson de las columnas predictivas con la variable ‘count’.

mean(p1_train$count)

[1] 191.5154

sd(p1_train$count)

[1] 182.139

3.1 Distribución

El grafico siguiente es un boxplot y un violin plot sobre la variable target ‘count’. Por el violin plot podemos ver que la mayoria de valores se encuentran cercanos al 0 pero por el boxplot vemos que los valores altos del dataset hacen que la media suba considerablemente. En el boxplot podemos apreciar también diversos outliers que sobresalen a partir del valor 550.

p1_train %>%
  ggplot( aes(x=count, y=count)) +
    geom_violin(width=1, color = "gray", alpha = 0.2, fill = 'green' ) +
    geom_boxplot(width=0.3, color="black", fill='blue', alpha=0.3, outlier.colour="red", outlier.shape=8,
             outlier.size=1, notch=TRUE) +
    theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
    ggtitle("Violin y box plot de los valores de la variable target ") +
    xlab("") + theme_classic() 

3.2 Histograma y densidad

Partiendo de la base que la variable target es continua, podemos apreciar que la frequencia de la misma sigue una distribución parecida a la exponencial, ya que cómo más se aleja del número 0 menos frequencia encontramos. Habrá que tener en cuenta entonces, que el modelo va a tender a predecir valores bajos en la mayoria de los casos.

p1_train %>%
ggplot( aes(x=count)) + 
 geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth = 20)+
 geom_density(alpha=.2, fill="#FF6666") + theme_classic()

4. Prediccion

4.1 Normalización y preparación de los datos

Hacemos un escalado sobre los datos para normalizar.

p1_train <- p1_train %>% mutate(across(where(is.numeric), scale))

4.2 Linear model

Entrenamos un modelo lineal

mod_lm <- lm(count~., p1_train) 

Calculamos el mse sobre el train

mean(mod_lm$residuals^2) #mse

[1] 19913.19

Calculamos el msle sobre el train

mean(log(mod_lm$residuals^2)) #mse

[1] 8.389741

Se puede apreciar que el mse y msle es muy grande, el modelo no funciona bien.

plot(p1_train$count)
abline(mod_lm)

Se puede ver cómo un modelo lineal es demasiado sencillo para entender la complejidad de los datos y senzillamente predice la media del resultado.

Podremos coger el error 20000 (que representa el mse entre la media y cada valor) para poder decidir en los modelos posteriores si estan aprendiendo.

4.3 glm

mod_glm <- glm(count~., p1_train, family = poisson())
mean(mod_glm$residuals^2)

[1] 0.7910352

Se ha intentado analizar la predicción pero no se encuentra sentido en la misma. El mse de 0.79 parece demasiado bajo para ser correcto y el predict sobre el train da valores count lejanos al ground truth.

4.4 Random forest

set.seed(1234) #definimos una random seed para obtener los mismos resultados consistentemente

Definimos el modelo. Durante la realización de la práctica se ha hecho un grid search y se definen los parametros que devuelven mejor resultado.

rf <- randomForest(count~., data = p1_train, mtry = 10, importance = TRUE, ntree = 50 )
rf

Call:
 randomForest(formula = count ~ ., data = p1_train, mtry = 10,      importance = TRUE, ntree = 50) 
               Type of random forest: regression
                     Number of trees: 50
No. of variables tried at each split: 9

          Mean of squared residuals: 2234.326
                    % Var explained: 93.26

Visualizamos que predictores son más importantes.

importance(rf)

              %IncMSE IncNodePurity
year       111.315166    22945095.4
hour       186.937499   154124280.2
season      40.132400     9248360.1
holiday      7.192508      699533.4
workingday 103.924054    19701681.3
weather     19.992407     3985389.6
temp        43.378082    30699917.8
humidity    27.163578     8386444.8
windspeed    6.956011     3337028.5

varImpPlot(rf)

Podemos apreciar cómo la variable ‘hour’ es con diferencia la que mejor ayuda a predecir el numero de bicicletas alquiladas.

Analisis del model performance (mse) sobre el conjunto de train

mseDF <- data.frame(pred = rf$predicted, gt = rf$y)
mseDF

Podemos ver las diferencias entre el count predecido y el ground truth.

mean((mseDF$pred-mseDF$gt)^2) #mse

[1] 2234.326

Calculamos RMSLE (error cuadratico medio logaritmico)

rmsle(mseDF$gt,mseDF$pred)

[1] 0.3695941

Obtenemos un mse y rmsle coherentes, el modelo funciona.

4.4.1 Random Forest 2

Volvemos a hacer random forest con variables más correlacionadas con la columna a predecir.

p1_trainImp <- p1_train %>% select('hour','year','workingday','season','humidity', 'count')

Entrenamos random forest

rf2 <- randomForest(count~., data = p1_trainImp, mtry = 10, importance = TRUE, ntree = 50 )
rf2

Call:
 randomForest(formula = count ~ ., data = p1_trainImp, mtry = 10,      importance = TRUE, ntree = 50) 
               Type of random forest: regression
                     Number of trees: 50
No. of variables tried at each split: 5

          Mean of squared residuals: 4122.14
                    % Var explained: 87.57

Comparación predicted vs ground truth.

mseDF2 <- data.frame(pred = rf2$predicted, gt = rf2$y)
mseDF2

Calculamos error

mean((mseDF2$pred-mseDF2$gt)^2) #mse

[1] 4122.14

Fucniona peor cogiendo solo las variables más importantes (mse 4122 > mse 2234). Nos quedaremos con la variable ’RF

5. Test

Realizaremos una predicción sobre el dataset de test con el mejor modelo, en nuestro caso el RF.

5.1 Realizar las mismas transformaciones en el test que en el train

p1_test <- p1_test %>% mutate(workingday = as.factor(workingday),
                                   holiday = as.factor(holiday),
                                   weather = as.factor(weather),
                                   season = as.factor(season)) %>% mutate(across(where(is.numeric), scale))

5.2 Realizar la prediccion

Se realiza la prediccion sobre los datos de test con el modelo RandomForest.

predsFin <- predict(rf, p1_test)
data.frame(prediction_test = predsFin)

3. Conclusiones

El algoritmo de regresion random forest es el que mejor predice el numero de bicicletas que se alquilaran un dia, con un mse de 2000 sobre el conjunto de train. Podemos ver analizando los resultados que el modelo no tiene mucho overfitting y generaliza bien, esto es asi en parte al haber creado un arbol no demasiado grande, evitando así la caída de valores en hojas.

Teniendo en cuenta que el mse entre la media y todos los valores del dataset era de 20.000 y con el mejor modelo hemos obtenido un mse de 2000, podemos concluir que nuestro modelo ha aprendido satisfactoriamente.

El rmsle definitivo es 0.3695941.

Cómo algoritmos de regresion se podrian intentar también el xgboost o el catboost pero perderiamos entendimiento de lo que esta haciendo el modelo, asi que se ha obviado para este trabajo.